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Although multilayer systems possess global inversion symmetry, some of the layers lack local inversion 
symmetry because no global inversion centers are present on such layers. Such locally non-centrosymmetric 
systems exhibit spatially modulated Rashba spin-orbit coupling. In this study, the superconductivity in mul- 
tilayer models exhibiting inhomogeneous Rashba spin-orbit coupling is investigated. We study the electronic 
structure, superconducting gap, and spin susceptibility in the superconducting state with mixed parity or- 
der parameters. We show the enhancement of the spin susceptibility by Rashba spin-orbit coupling and 
interpret it on the basis of the crossover from a centrosymmetric superconductor to a non-centrosymmetric 
superconductor. It is also shown that the spin susceptibility is determined by the phase difference of the 
order parameter between layers and is nearly independent of the parity mixing of order parameters. An 
intuitive understanding is given on the basis of the analytic expression of superconducting order parameters 
in the band basis. The results indicate that not only a broken global inversion symmetry but also a broken 
local inversion symmetry leads to unique properties of superconductivity. We discuss the superconductivity 
in artificial superlattices involving CeCoIus and multilayer high-Tc cuprates. 
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1. Introduction 

The discovery of superconductivity in the heavy- 
fermion compound vifithout inversion symmetry, 
CePtaSi,^^ triggered extensive studies of non- 
centrosymmetric superconductivity.^-' Subsequently, 
several new non-centrosymmetric superconductors with 
unique properties have been identified among heavy- 
fermion materials and in other classes of materials. 
The research field has even been effectively extended 
to non-centrosymmetric superfluids in cold Fermi gases 
18, 19) jjj^ which antisymmetric spin-orbit coupling is arti- 
ficially induced. Many interesting properties, such as the 
parity mixing of order parameters,^*^-* the magnetoelec- 
tric effect,^^"^**^ anisotropic spin susceptibility^^' ^^"^^^ 
accompanied by an anomalous paramagnetic depairing 
effect, ^'^^ helical superconducting phases in a magnetic 
field, '^^"^^^ and topological superfluid phases, have 
been proposed and studied in various contexts. 

The term non-centrosymmetric superconductivity is 
commonly used for systems that have no inversion cen- 
ter. It is, however, interesting to note that features fa- 
miliar in non-centrosymmetric superconductors can also 
be relevant for systems with an inversion center, but 
specific forms of a local violation of inversion symme- 
try. A clear example is a multilayer system in which the 
local mirror symmetry is broken, as shown in Fig. I. 
Naturally, the question arises whether superconductiv- 
ity displays any exotic property in such a locally non- 
centrosymmetric system. One of the authors has earlier 
investigated this issue by studying the spin triplet su- 
perconducting state with random Rashba spin-orbit cou- 
pling arising from stacking faults. In this study, we 
focus on regular multilayer systems and elucidate the 



basic properties of such locally non-centrosymmetric su- 
perconductors. We show that, indeed, some properties of 
the superconducting phase can be strongly affected by 
a broken local inversion symmetry in multilayer systems 
that exhibit small interlayer coupling. 
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Fig. 1. (Color online) Schematic view of trilayer system with 
focus on local inversion symmetry. The inversion center is present 
in the inner layer; it is absent in the outer layers. This crystal 
structure is regarded as a locally non-centrosymmetric system. 



Our study is motivated by the observation of super- 
conductivity in artificial superlattices consisting of the 
heavy-fermion superconductor CeCoIns and the conven- 
tional metal YbCoIns.^^'^^^ Although it is expected that 
the superconductivity occurs in sufficiently thick multi- 
layer structures of CeCoIn5, it is surprising that it even 
prevails down to superlattice incorporating stacks of only 
three well-separated CeCoIns layers. 

The superconductivity in bulk CeCoIns has attracted 
considerable interest for its unique magnetic features, 
such as the paramagnetic depairing effect, ^"^^ the pos- 
sible realization of the Fulde-Ferrell-Larkin-Ovchinnikov 
phase at high magnetic fields,^**' an unconventional 
magnetic order, '*^''*^^ and field-induced quantum criti- 
cality.'*^"^^^ These features make it even more attractive 
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CeCoIno. As such, it represents an ideal system for study- 
ing the effects of local inversion symmetry breaking, and 
more so because many of the striking properties of non- 
ccntrosymmetric superconductors may result from the 
additional presence of strong magnetic correlations. 

A better class of systems of similar character is that 
of multilayer high-Tc cuprates.^'^"^^^ Although cuprate 
superconductors have been intensively investigated since 
their discovery in 1986, the role of spatially modulated 
Rashba spin-orbit coupling has not received much atten- 
tion thus far. Since the magnetic properties of multilayer 
cuprates are investigated by nuclear magnetic resonance 
(NMR) measurement, it is also possible to study multi- 
layer cuprates focusing on local non-centrosymmetricity. 
For instance, the spin susceptibility in the superconduct- 
ing state can be measured with a layer resolution us- 
ing the Knight shift of NMR. Following a brief report 
on the spin susceptibility in locally non-centrosymmetric 
systems, we study here the electronic structure, su- 
perconducting gap, and spin susceptibility of multilayer 
superconductors in detail. 

The paper is organized as follows. In §2, we intro- 
duce the model Hamiltonian for multilayer systems in- 
cluding spatially modulated Rashba spin-orbit coupling 
and parity-mixed superconducting order parameters. In 
§3, we show the numerical results of the spin suscep- 
tibility in the bi- and tri-layer systems to demonstrate 
the crossover from a conventional superconductor to a 
non-centrosymmetric superconductor. The numerical re- 
sults arc discussed on the basis of analytic expressions 
for the electronic structure (§4) and superconducting gap 
(§5) by decomposing the spin susceptibility into the Pauli 
and Van-Vlecks parts (§6). The relationship between the 
symmetry of order parameters and the spin siisccptibility 
is clarified in §7, where a nontrivial role of the interlayer 
phase difference of order parameters is discussed. The 
numerical results of the spin susceptibility in more than 
three layers are shown in §8. A brief summary and dis- 
cussion of the superconductivity in CeCoIns are given in 
§9. 

2. Model 

First, we introduce a model Hamiltonian for two- 
dimensional multilayer superconductors with spatially 
modulated Rashba spin-orbit coupling as 

H = + Hsoc + Hpaii + H±, (1) 

^b=Yl ^(^) 4>sm(^ksm, (2) 

HsOC= ^ am9{k) ■ a-ss'clsrn^ks'm, (3) 

k.s,s' ,m 

Hp^ir = \ Yl l^ss'm{k)cl^^cl^,,^ + h.C.], (4) 

k,s.s' ,m 

H±=t± 4smCfesm', (5) 

where Cksm {Cksm) annihilation (creation) operator 
for an electron with spin s on layer m, and cTss' is the 



The first term, ffb, describes the dispersion relation 
without spin-orbit coupling or interlayer coupling. We 
consider a square lattice and assume a tight-binding 
model, i.e., e{k) = — 2t{cos kx+ cos ky) — fi. The {x,y,z)- 
axes correspond to the (a, b, c)-axes of the tetragonal 
crystal structure. We choose the unit of energy as t = 1 
and fix the chemical potential = — 1 , which leads to the 
electron density per site being approximately 0.63. The 
following results are nearly independent of the dispersion 
relation and electron density. 

The second term, Hsoc, describes the Rashba spin- 
orbit coupling arising from the lack of local inversion 
symmetry. This term preserves time reversal symmetry, 
if the g-vector is odd in fe, i.e., g{—k) = —g{k). The 
coupling constants am should have opposite signs be- 
tween layers above and below the inversion center so 
as to conserve global inversion symmetry. For instance, 
(ai, a2, as) = {a, 0, —a) for trilayers. For g{k), the mul- 
tilayer systems should have a Rashba-type g-vector be- 
cause the mirror symmetry is broken for outer layers (see 
Fig. 1). Although the detailed momentum dependence of 
g{k) is determined by electronic structures, we assume 
here the simple form g{k) — (— sin fcj,, sin fc^;, 0). 

The third term, TJpair, introduces intralayer Cooper 
pairing via an off-diagonal mean field. We ignore inter- 
layer pairing as we consider the layers to be weakly cou- 
pled. Since we take into account the spatially modulated 
Rashba spin-orbit coupling arising from a broken local 
inversion symmetry, the order parameter Ass'mO^) 
volves both spin singlet and triplet components, 

A — f -d^m{k) + idymik) ipmik) + dzmik) \ 

(6) 

where tpm{k) and dm{k) are the corresponding scalar 
and vector order parameters for the spin singlet and 
triplet pairings on the layer m, respectively. For our 
discussion, we avoid the use of a microscopic model 
based on a pairing mechanism. Rather, we introduce 
an order parameter on phenomenological grounds, which 
we assume to have s-wave symmetry for the singlet 
pairing and a p-wave symmetry for the triplet pairing. 
On symmetry grounds, we choose ijjm{k) = ipm and 
dmik) = dmg{k) = dm{-sm.ky,s.m.kx,Q). We take |V'm|, 
\dm\ < 0.01 to be sufficiently small to satisfy the condi- 
tion |Ass'm(fc)| ^ jam I ^ ep (ep is the Fermi energy), as 
realized in most non-centrosymmetric superconductors. 
To minimize the interlayer coupling energy, the dominant 
order parameter component maintains the same sign over 
all layers, whereas the other (subdominant) component 
adjusts the sign to that of spin-orbit coupling (a^)- 

The fourth term, describes the interlayer hopping 
of electrons between nearest-neighbor layers. Since we 
consider a quasi-two-dimensional system, we assume that 
the interlayer hopping t±_ is smaller than the intralayer 
hopping t. 
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3. Numerical Results of Spin Susceptibility 

3.1 Spin susceptibility for fields along c-axis 

We now determine the spin susceptibility of multilayer 
superconductors for spatially modulated Rashba spin- 
orbit coupling, with a magnetic field applied along the 
c-axis. The spin susceptibility x = limff-i-o(-^^s)/-ff is ob- 
tained by calculating the magnetization (Ms) in the field 
H. The Zeeman coupling term is introduced as 

Hz = -^^ ^ H ■ (Tss'cl^,n(^ks'nu (7) 

where we assume g = 2 and /ie is the Bohr magneton. 
First, the Hamiltonian is diagonalized by introducing the 
unitary transformation (7^ = f^C/^(fe) in the Nambu 
space of M layers, where quasi-particle operators form 
a 4M-dimensional vector, 

' ■ • ■ ' Cfc-fM' "^fciA/' C-fetM , C-klAi), (8) 

and are analogous to the Bogoliubov quasi-particle 
operators = (tIu 7fe2' ■ • • > 7fc4M)- The diagonalized 
Hamiltonian is written as 

^ 4M 
k i=l 

where Ei{k) is the quasi-particle energy. The magnetiza- 
tion is given as 

4M 

{Ms} - ^EE[^'(^)]^^/(^»w)' (10) 

k i=\ 

where /(-E) is the Fermi-Dirac distribution function. The 
matrix representation of a spin operator is defined on the 
t\ basis as 

= ij\k)t^t]{k), (11) 

with E'' as the /x-component of the spin operator in the 
4M-dimensional space. 

We now focus on the spin susceptibility of bilayer 
(M = 2) and trilayer (M = 3) systems at zero tem- 
perature. More than three layers will be discussed in §8. 
The corresponding coupling constants of Rashba spin- 
orbit coupling are described as (ai,a2) = (a, —a) for 
bilayers and (ai,a2,a3) — (a,0, — a) for trilayers. 

We compare two cases: (A) the spin triplet channel is 
dominant |d| > 1-01 and (B) the spin singlet channel is 
dominant > \d\. The order parameters are assumed 
as follows: 

(A) (V'l, "02) — {ip, —ip) and (^1,^2) — {d, d) for bilayers. 
(V'i,i/'2,i/'3) = (■0,0, -■0) and {di,d2,dz) = {d,d,d) 
for trilayers. 

(B) (01,02) = {ip,i^) and (^1,^2) = {dj—d) for bilayers. 
(01, "02, "03) = (0',V','0) and {di,d2,d3) = {d,0,-d) 
for trilayers. 

In case (A) , the spin susceptibility remains unaffected by 
the superconductivity Xs = Xn because the spin triplet 
component of the type dm{k) cx g{k) _L z is an equal- 
spin pairing state with Cooper pair spins along the c-axis. 



possible without pair breaking. This feature is essentially 
independent of spin-orbit coupling, interlayer hopping, 
and the number of layers, as can be seen in Figs. 2 and 
3 for both bi- and tri-layer systems, respectively. We will 
provide more rigorous arguments for this kind of behav- 
ior in §7. 
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Fig. 2. (Color online) Spin susceptibility along c-axis for bilayer 
system. The interlayer coupling t± is shown in the figure for the 
dominantly spin singlet pairing state (dashed, dotted, and dash- 
dotted lines). The spin susceptibility is independent of interlayer 
coupling when the spin triplet channel is dominant (solid line) . We 
assume that ip = 0.01 and d = in the former, and i/, = and 
d = 0.01 in the latter. The spin susceptibility is nearly independent 
of the subdominant order parameter, as will be shown in Fig. 5. 
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Fig. 3. (Color online) Spin susceptibility along c-axis for trilayer 
system. The interlayer coupling t± is shown in the figure for the 
dominantly spin singlet pairing state (dashed, dotted, and dash- 
dotted lines). The spin susceptibility is independent of interlayer 
coupling when the spin triplet channel is dominant (solid line). The 
spin susceptibility is nearly independent of the subdominant order 
parameter as in the bilayer system (Fig. 5). 

More interesting is case (B), as spin singlet pairing 
leads to a complete suppression of the spin susceptibility 
at T = in a conventional superconductor. Indeed for a 
vanishing spin-orbit coupling (a = 0), we find Xs = ir- 
respective of t± . As soon as spin-orbit coupling is turned 
on, however, the spin susceptibility gradually recovers 
and approaches a constant value for large a: Xs — ^ Xn for 
the bilayer and Xs ^ 2xn/3 for the trilayer. The mecha- 
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by Rashba-type spin-orbit coupling (see §4). Note that 
the layers behave as being nearly decoupled for large a, 
and then the spin susceptibility along c-axis is recovered 
for layers with non- vanishing am- Consequently, in the 
bilayer system, all layers are involved, giving rise to a full 
recovery of Xs for large a (analogous to the uniformly 
non-centrosymmetric superconductor^®^ ) , whereas in the 
trilayer system, only two of the three layers can con- 
tribute to what yielding a correspondingly reduced lim- 
iting value of Xs — 2xn/3. Figure 4 corroborates this pic- 
ture by considering the contributions of the different lay- 
ers. Indeed, in a large-a regime, the outer layers m = 1, 3 
carrying spin-orbit coupling saturate at Xs Xn/3, while 
the center layer m = 2 completely suppresses the spin 
susceptibility. Remarkably, at small a {< t±), Xs behaves 
in the same way for all layers and surprisingly leads to a 
nonmonotonic a-dependence for the center layer. 




Fig. 5. (Color online) Zero-temperature spin susceptibility along 
c-axis for various \ip\ values in bilayer system. We assume the am- 
plitude of the spin triplet component to be \d\ = 0.01— We show 
the results for interlayer coupling t± = 0.1 and spin-orbit coupling 
a = (crosses), a = 0.05 (triangles), and a = 0.2 (diamonds). 




a 



Fig. 4. (Color online) Contribution of each layer to spin suscep- 
tibility along c-axis in trilayer system. We assume the spin singlet 
state {i> = 0.01 and d = 0) with = 0.1 (solid line) and tj_ = 0.2 
(dashed line). 



the interlayer phase differences of order parameters (§7). 
Since we assume the zero-phase difference for the domi- 
nant order parameter and the 7r-phase difference for the 
subdominant one, the spin susceptibility is determined 
by the dominant order parameter and is negligibly af- 
fected by the subdominant component. The details will 
be given in §7. 

3.2 Spin susceptibility along ab-axis 

Within our model, we find that the spin susceptibility 
along the a&-axis is always half of that along the c-axis, 
independent of the strength of a and t± as well as of the 
number of layers. For numerical evidence. Fig. 6 shows 
the spin susceptibility along the a6-axis in the trilayer 
system; it is one-half of the results in Fig. 3. We con- 
firmed these behaviors for another number of layers. 



The numerical data in Figs. 2 and 3 show that in- 
terlayer hopping is in competition with spatially modu- 
lated Rashba spin-orbit coupling, such that a larger t± 
yields a higher effective acfr ^ t±_ for the crossover from 
the behavior of a conventional superconductor to that of 
a non-centrosymmetric superconductor. For aeS 3> tj_, 
layers are almost decoupled from each other (see §4), 
and then the multilayer superconductor is regarded as 
a set of (non-)ccntrosymmetric superconductors. This 
crossover is best observed in the peak of Xs a-t approx- 
imately Qfoff ^ t± for the center layer of the trilayer 
system (Fig. 4). Thus, modifying t±, e.g., by applying 
uniaxial stress along the c-axis, can affect the magnetic 
response for c-axis fields in case (B). No such effect is 
expected in case (A). 

Surprisingly, these results are almost independent of 
the subdominant order parameter. Figure 5 shows the 
spin susceptibility of bilayers for various values while 
keeping the summation \tfj\ + \d\ = 0.01. We see a 
nearly constant spin susceptibility except for the jump 
at IV'I = \d\ = 0.005. This jump arises from a discontinu- 
ous change of the order parameter, since we assume case 
(A) for I'f/'l < \d\ and case (B) for \ip\ > \d\. We will show 
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Fig. 6. (Color online) Spin susceptibility along a-axis for trilayer 
system. Parameters are the same as those in Fig. 3. 



4. Electronic Structure in Normal State 

To elucidate the crossover from the behavior of 
a centrosymmetric superconductor to that of a non- 
centrosymmetric superconductor shown in the numerical 
calculation (Figs. 2 and 3), we give an analytic expression 
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gap in §4 and §5, respectively. An intuitive understand- 
ing is given in §6 by decomposing the spin susceptibility 
into the Pauli and Van-Vleck parts. 

4-1 Bilayer system 

First, we study the single-particle wave function in the 
normal state of bilayer systems. The schematic figure for 
the role of Rashba spin-orbit coupling and interlayer cou- 
pling is shown in Fig. 7. When the two layers are decou- 
pled at t± ~ 0, each layer has Fermi surfaces, as shown 
in Fig. 7. The structures of Fermi surfaces are the same 
for layers 1 and 2, but the spin orientations are opposite, 
because the Rashba spin-orbit coupling am has the op- 
posite sign ai = —ai. When the interlayer coupling 
is switched on, the Fermi surfaces are coupled, as shown 
by dashed arrows in Fig. 7. 




m=l 



m=2 



Fig. 7. Schematic figure of Fermi surfaces in bilayer system witfi- 
out interlayer coupling tx_ = 0. In this case, layers are decoupled, 
and then the Fermi surfaces of layers 1 and 2 are independent 
of each other. The Fermi surfaces are split by Rashba spin-orbit 
coupling in both layers 1 and 2, but the spin orientations are op- 
posite because Rashba spin-orbit coupling am, has the opposite 
sign (ai = — 02). When interlayer coupling is tuned on, the 
electronic states are coupled, as shown by dashed arrows. 



We show the wave function of quasi-particlcs with the 
momentum k — (0, fcy) as an example. According to the 
Fermi surfaces shown in Fig. 7, a simple expression is ob- 
tained by choosing the spin quantization axis along the 
a;-direction. In the following presentation, the wave func- 
tions I —i')m and I describe the electron states on the 
layer m with right- and left-pointing spins, respectively. 
In the normal state with Ass'm{k) = 0, the Hamiltonian 
[eq. (1)] is block-diagonalized in this representation. The 
block of the Hamiltonian is obtained as 



e{k)+a'{k) 



and 



e{k)-a'{k) 



e{k)-a'{k) 



e{k)+a'{k) 



for 



for 



(12) 



, (13) 



respectively. We denote the magnitude of spin-orbit cou- 
pling as 



with twofold degeneracy: 
^^^\k)^e{k) 



E 



a'iky+tl, 



Ei^^ (k) = e{k) ~ ^a'ik)^+t] 



The wave function is obtained as 



v^-Aikr\ - ^)i+A(fc) i 
A{k)\^), + ^i-A{kr\ 



for E ^ E[^\k) and as 



- A(fc)|^)i + Vl-A(fc)2| 

v^-Mkr\ ^)i-A(fc)i^)2 

(2) 

for E = £'2 (fe). We here defined 
A{k) 



^tl+(a'{k) + ^a'{kY+tiy 



(15) 
(16) 



(17) 



(18) 



(19) 



a' [k) = a\g{k)\ — aJ sm fc^ + sin ky. 



(14) 



We obtain a similar single-particle wave function for the 
other momentum fcx 7^ by choosing a spin quantization 
axis parallel to g{k). 

The twofold degeneracy in the above two energy bands 
originates from the time-reversal and global inversion 
symmetry. Thus, the structure of the energy bands 
seems to be the same as that in conventional metals. 
However, Rashba spin-orbit coupling affects the single- 
particle wave function in a unique way. Note that A{k) in 
eq. (19) is determined by the competition between inter- 
layer coupling t±_ and spin-orbit coupling a. When spin- 
orbit coupling is absent, i.e., a — eq. (19) is reduced 
to A{k) = 1/^/2. Then, two energy bands are regarded 
as bonding and antibonding bands with spin degeneracy, 
namely, {l/^/2){\ ^)i + | ^)2) and (l/^/2){\ + 
)2) for E = E^^\k), and (1/V2)(| ^)i - | ^)2) and 
(1/V2)(| ^)i - I <-)2) for E = E''^\k), respectively. 
This is the conventional band structure of bilayer sys- 
tems. On the other hand, we obtain A{k) = in the 
limit of large spin-orbit coupling a ^ t±, and then the 
twofold degeneracy is given by | — t^) 1 and 1^)2, and by 
I and I — 7>)2- Thus, the quasi-particle is localized on 
each layer and degenerates with a quasi-particle of oppo- 
site spin on the other layer. In this case, the splitting of 
the two energy bands is induced by the spin splitting in 
each layer, although the twofold degeneracy is protected 
by the global inversion symmetry. With increasing spin- 
orbit coupling, the two energy bands change their charac- 
ter from the bonding and antibonding states to the spin 
split states on each layer. This crossover leads to the a 
dependence of spin susceptibility in Fig. 2. 

4-.2 Trilayer system 

Next, we study the trilayer systems. The single-particle 
wave function is obtained in the same way as that in the 
bilayer system. Figure 8 shows the schematic figure of 
Fermi surfaces in the decoupling limit t^ — Q. For the 
quasi-particles with k = (0, fcy), the Hamiltonian without 



Diagonalizing the 2x2 matrix, we obtain eigenvalues 
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m=l 



m=2 



m=3 



Fig. 8. Schematic figure of Fermi surfaces in trilayer system 
without interlayer coupling t±_ = 0. Fermi surfaces in layers 1 and 3 
are split by Rashba spin-orbit coupling; however, they degenerate 
in layer 2 because of Q2 = 0. When interlayer coupling is turned 
on, the electronic states are coupled, as shown by dashed arrows. 



pairing field Ass'm{k) = is block-diagonalized as 
'e{k) + a'{k) tj_ \ /| ^)i 



t± s{k) t±^ I for 
tj. e{k)-a'{k) 




and as 



'e(fc) - a'ik) 

t± e{k) I for I I ^) 2 I .(21) 

e{k)+a'{k)) VI ^)3, 

Diagonalizing the 3x3 matrix, we obtain tliree eigen- 
values: 



EY\k)^e{k) + \la'{kY + 2t 
E^^\k)^e{k), 



(22) 
(23) 



Ef\k)^e{k)-^a'{kY + 2tl, (24) 

with twofold degeneracy. The wave function is described 
as 



SI Hi + C\ ^2 + ^/l-B^-C^ ^3 
for E = Ef\k), 

-C\ Hi + Vl-2C2| H2 + C| Hs 



(25) 



(26) 



(27) 



C\ ^>i+Vl-2C2| ^>2-C| Ha ' 

for E = E!f\k), and 

B\ Hi-g| + Vl-B^-C^\ Hs 
Vi-B2_c2| Hi - c| H2 + s| Hs 

for E — E^ (fc)- We omitted the index fc in eqs. (25)- 
(27), and defined B{k) and C{k) as 

^2 

= ^ , , 

(a'ik) + ^a'ikf + 2t^)^a'(fc)2 + 2t^ 

C(fc): * 



, (28) 
(29) 



Va'(fc)2 + 2i] 
respectively. 

In the absence of spin-orbit coupling (a = 0), we ob- 
tain B{k) = 1/2 and C{k) = l/\/2; then the energy 
bands with E = E\ (k) and E = E^ (fc) have an even 



parity with respect to the center layer, while the band 
E = E^\k) has the odd parity. When the spin-orbit 
coupling is turned on, the parity is mixed. In the limit 
of large spin-orbit coupling a ^ t±, the quasi-particle 
is localized on each layer, as shown in Fig. 8. The dou- 
bly degenerate wave functions are | — ^)i and | -s— )3 for 
E = Ef\k), I ->)2 and | ^)2 for E = E^^^\k), and 
I Hs and I Hi for E = E^f\k). 

5. Order Parameter in Band Basis 

We turn to a superconducting order parameter in the 
band basis. For this purpose, the Hamiltonian is trans- 
formed by the unitary matrix T(fe), which diagonalizes 
the Hamiltonian without a pairing field {Ass'm{k) = 0), 



H 



\ CiH{k)Ck = \Y. ZlH'{k)Zk. (30) 
fc fc 



The basis of zj^ is chosen to be the eigenstates in 
eqs. (17) and (18) for bilayers. Then, we obtain the fol- 
lowing matrix representation of the Hamiltonian; 

H'{k) = f^k)H{k)f{k) 
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(31) 

where the index k is omitted for simplicity. 

Since we consider a small superconducting gap |A| ^ 
\a\, as realized in most (locally) non-centrosymmetric su- 
perconductors, the role of interband Cooper pairing is 
ignored. Then, the superconducting order parameter of 
each band is described by Ai(fc) because the 2x2 ma- 
trix for intraband Cooper pairing is diagonalized in this 
representation. Thus, the superconductivity in bilayers 
is regarded as an equal pseudo-spin pairing state on the 
basis of eqs. (17) and (18). 

The superconducting order parameter in each band is 
obtained as follows. There are two cases with respect to 
the interlayer phase difference, as discussed in §2. In case 
(A), the spin triplet order parameter (d- vector) has the 
same sign in both layers 1 and 2, while the spin sin- 
glet order parameter has opposite signs in the two lay- 
ers, namely, {ipi,ip2) = ii',—^') and {di,d2) — {d,d). 
In case (B), the spin singlet order parameter has the 
same sign in layers 1 and 2, but the d-vector has oppo- 
site signs in the two layers, namely, (V'l, "02) — i') and 
(di,d2) = {d,-d). 

Case (A) leads to the order parameter in the band 
basis for k — (0, fcy): 



Ai(fc) = -fc_ 



Va2|g(fc)|2+4 



i/j + d 



(32) 
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A2(fe) = Ai(fe), 

A3(fe) = -fc_ 



A4(fc)=A3(fe), 



(33) 

(34) 



(35) 



where k± = smky ± isinfcj;. On the other hand, we 
obtain the fohowing results in case (B): 



Ai(fe) = 



\9{k)\ 
A2(fe) = -Ai(fc), 

A3(fe) ^ 



l9(fe)l 



\/«%(fc)P 



a\g{k) 



A4(fc)--A3(fe). 



(36) 

(37) 
, (38) 
(39) 



Thus, the relative sign of Ai(fc) and A2(A;) and that of 
A3(fe) and A4(fe) are different in cases (A) and (B) and 
thus also yield different spin susceptibilities, as shown in 
Fig. 2 and discussed in §7. We do not show the super- 
conducting order parameter in the trilayer systems, but 
it is straightforward to extend our analytical calculation 
to trilaycrs. 

6. Spin Susceptibility 

6.1 Pauli and Van-Vleck contributions 

We provide here the theoretical basis of our numerical 
results of spin susceptibility in §3. For this purpose, wc 
decompose the spin susceptibility in the normal state into 
the Pauli and Van-Vleck parts. With the use of the linear 
response theory, the transverse component of uniform 
spin susceptibility is obtained in the normal state as 



f{E,{k))-f{E^{k + q)) 
E^{k + q)-Er,{k) ' 

which is decomposed into 



(40) 



Xn 



The Pauli spin susceptibility arises from the in- 
traband contributions, while the interband polarization 
gives rise to the Van-Vleck spin susceptibility ■ They 

can be expressed as 



(42) 



X 



J2 Y.(^\S+\u){u\S-\v)6{Er,{k)), 

E^=E„ k 



Er,i^E„ k 



E^k) - Enik) 



respectively. In multilayer systems, the Fermi surfaces 
are split by both spin-orbit coupling and interlayer cou- 
pling, whereby the former gives rise to the Van-Vleck 
susceptibility. 

For the discussion of the spin susceptibility on the basis 

of the single-particle wave functions given in §4, it should 

Vip nnf.pH that thp snin nni^nti^i^tinn piyis is HirprtpH alnncr 



the a6-planc in the entire Brillouin zone. Since the spin 
quantization axis is perpendicular to the c-axis, the spin 
susceptibility along the c-axis is given by the transverse 
spin component that can be calculated using cq. (40). 

The Van-Vleck spin susceptibility is negligibly affected 
by the superconducting gap when the band splitting due 
to the spin-orbit coupling and interlayer coupling is much 
larger than the superconducting gap. Thus, the Van- 
Vleck contribution provides a lower limit of spin suscep- 
tibility in the superconducting state. On the other hand, 
the Pauli spin susceptibility may be suppressed by the 
superconductivity and depends strongly on the symme- 
try of the order parameter. We will show that the Pauli 
spin susceptibility of multilayer superconductors is de- 
termined by the phase difference of the order parameter, 
which has been discussed for bilayers in §5, irrespective 
of the ratio of the spin singlet and the triplet order pa- 
rameters. The Pauli spin susceptibility is completely sup- 
pressed in case (B), whereas it remains unaffected in case 
(A), as will be shown in §7. 

We assumed case (A) [case (B)] in the dominant spin 
triplet (singlet) pairing state for the numerical calcula- 
tion in §3. Thus, the spin susceptibility in the dominantly 
spin singlet superconducting state studied in §3 coincides 
with the Van-Vlcck part in the normal state, and there- 
fore Xs/Xn = X"' I (X^ + X^), but Xs/Xn = 1 in the domi- 
nantly spin triplet superconducting state. We discuss the 
Pauli and Van-Vleck spin susceptibilities of bi- and tri- 
laycrs below. 

6.2 Bilayer system 

With the use of the single-particle wave functions 
eqs. (17) and (18) in the bilayer system, we obtain the 
normal state Pauli and Van-Vleck spin susceptibilities at 
T = as 



X 



J2 (1 - Aikf) [6iE['^ (fc)) + <5(ijf ) (fc))] 



p{0){4A{kf{l-A{kf) )fs, 



(41) and 



^2(l-2A(fc) 



e{Ei'\k)) 



&iEi'\k)) 



Ei'\k)-El'>{k) 



(2), 



p(0)((l 



2A(fc)2)')ps, 



(44) 
(45) 

(46) 
(47) 



respectively. Note that the Fermi-Dirac distribution func- 
tion f{E) is reduced to the step function Q{E) at T = 0. 
The average on the Fermi surface is denoted as (• • • )fs 
and p(0) is the density of state at the Fermi level. 

In the absence of spin-orbit coupling, the Van-Vlcck 



(43) P^'"* vanishes and the Pauli part is obtained as x 



because of A{k) = 1/^/2. On the other hand, the Pauli 
part vanishes for t± = such that x^ = X^~ because of 
A{k) = 0. 

We show the numerical results of x^/(x^ + X^) ob- 
tained using eqs. (44) and (46) for various spin-orbit cou- 
plings a in Fig. 9. We see that the Van-Vleck part of the 
spin susceptibility in the normal state, X^/(x^ + X^); 
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Fig. 9. (Color online) Van-Vleck part of spin susceptibility along Fig. 10. (Color online) Van-Vleck part of spin susceptibility 
c-axis in normal state of bilayer system, x^/ (x^ +X^) (solid lines). along c-axis in normal state of trilayer system, x^/ (x^ + X^) (solid 



The dashed lines show the spin susceptibility in the superconduct- 
ing state, Xs/Xn> for ip = 0.01 and d = 0, which is shown in Fig. 2. 
We see that these quantities coincide with each other. The param- 
eter is chosen as t± = 0.1 or 0.2. 



lines). The dashed lines show the spin susceptibility in the super- 
conducting state, Xs/Xn> for ip = 0.01 and d = 0, which is shown 
in Fig. 3. We see that these quantities coincide with each other, as 
in the bilayer systems. The parameter is chosen as t± =0.1 or 0.2. 



bility, Xs/Xn, in the superconducting state with a dom- 
inant spin singlet order parameter. Thus, we obtain an 
estimate of the spin susceptibihty in the superconduct- 
ing state through the Van-Vleck spin susceptibility in the 
normal state. 

6.3 Trilayer system 

With the use of the single-particle wave functions 
eqs. (25)-(27) in the trilayers, the normal state Pauli spin 
susceptibility at T = is obtained as 



X 



+ {l-4C^yS{Ei'^^)\ (48) 



2p(0) 



+ 



3 
P(0) 



FS 



FS, 



(49) 



X 



while the Van-Vleck spin susceptibility is obtained as 



k 
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(50) 



(2C2(Vl - 52 - C2 + - 2C2 - B) 



;fs 



+ ^(2(1 -2C^)^) 



FS, 



(51) 



where we omitted the index k in E^ (fc), B{k), and 
C{k) for simplicity. 

The Van-Vleck spin susceptibility vanishes in the ab- 
sence of the spin-orbit coupling a = 0, and x^ ~ 
2p(0)/3 = 2xn/3 in the limit of large spin-orbit coupling 



a 3> t±, as expected from the discussions in §3 and §6.1. 
Figure 10 shows that the Van-Vleck part of the spin sus- 
ceptibility in the normal state x^/(x^ + X^) given by 
eqs. (48) and (50) coincides with the spin susceptibility 
in the superconducting state Xs/Xn, which is shown in 
Fig. 3. Figure 10 shows that the Pauli spin susceptibil- 
ity is completely suppressed in the dominant spin singlet 
pairing state; therefore, the spin susceptibility of trilay- 
ers at r = is given by the Van-Vleck part, as in the 
bilayer system. 

7. Relation to Symmetry of Superconductivity 

In this section, we examine the relation between the 
order parameter and spin susceptibility in our multilayer 
superconductors. As we discussed in §6, the Van-Vleck 
part of the spin susceptibility is not affected by super- 
conductivity. On the other hand, the Pauli spin suscep- 
tibility has a clear dependence on the superconducting 
order parameter. Indeed, as we will show later, the Pauli 
spin susceptibility is determined by the phase difference 
of the order parameter between layers, and is essentially 
independent of the ratio of parity mixing, 

For the bilayer system, we consider cases (A) and (B), 
discussed in §5, to optimize the interlayer Josephson cou- 
pling energy. 

To understand the above-described aspect of the order 
parameter structure, it is important to view the situation 
in the band basis, as in eqs. (32)-(39). For fe = (0, ky), in 
the band basis, the "d- vector" of the pseudo-spin points 
along the y-axis in case (A) because Ai(fc) = A2(fc) 
and A3(fc) = A4(fe). On the other hand, in case (B) the 
"d- vector" is oriented along the .x-axis, whereby we find 
Ai(fc) = -A2(fc) and A3(fc) = -A4(fc). Since the spin 
quantization axis is along the a-axis of tetragonal crys- 
tals for fc = (0, fcy), the crystal c-axis corresponds to the 
y-axis with the z-axis being the spin quantization axis. 
Thus, the "d- vector" of a pseudo-spin is parallel to the 
c-axis in case (A) and perpendicular to the c-axis in case 
(B). This is true not only for k = (0, fcy) but also for all 
momenta in the Brillouin zone. One may mistakenly con- 



Ti — l; . 
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c-axis is suppressed (unchanged) in case (A) [case (B)] 
similarly to a spin triplet superconductor with a d- vector 
parallel (perpendicular) to the field. However, the oppo- 
site relation is actually obtained for the " d- vector" in the 
band basis, because the spin of each band is inversely po- 
larized for k and ~k. Therefore, taking into account the 
opposite sign of the spin quantization axis, it is shown 
that the Pauli spin susceptibility along the c-axis remains 
unchanged by the superconductivity in case (A), while it 
is completely suppressed in case (B). 

We conclude here that the spin susceptibility is deter- 
mined by the interlayer phase difference of the order pa- 
rameters and independent of their amplitude. This sur- 
prising result has already been anticipated from Fig. 5. 
The spin susceptibility is nearly constant for ip, although 
it shows a jump at \^\ = \d\, following the discontinuity 
of the phase difference. Note that we assumed case (A) 
for IV^I < \d\ and case (B) for lip] > \d\. 

The above finding is confirmed by calculating the spin 
susceptibility while assuming case (A) or (B) regardless 
of the amplitude of singlet and triplet order parameters. 
Figure 11(a) shows the result in case (A), while Fig. 11(b) 
shows that in case (B). It is shown that the spin suscep- 
tibility at T = is nearly independent of the ratio of 
parity mixing from = to = oo. Thus, the 

spin susceptibility remains unchanged when the domi- 
nant order parameter is a spin singlet component with 
a TT-phase difference. On the other hand, the Pauli spin 
susceptibility is completely suppressed by the spin triplet 
order parameter with a 7r-phase difference between lay- 
ers. In other words, the spin singlet order parameter with 
a TT-phase difference plays the role of a spin triplet order 
parameter with a zero-phase difference, and vice versa. 

We obtain the same conclusion for the trilayer sys- 
tems. The spin susceptibility at T = is determined by 
the phase difference and is independent of the amplitude 
of the order parameters. We confirmed that the spin sus- 
ceptibility of trilayer superconductors is nearly constant 
for IV'I/MI, as in Fig. 11. 

8. More than Three Layers 

Finally, we show the numerical results of the spin sus- 
ceptibility for systems with 4,5,6, and 7 layers. Su- 
perconductivity has been observed on the artificial su- 
perlattice of CeCoIus with 3, 4, 5, 7, and 9 layers. We 
consider the dominantly spin singlet pairing state with 
= 0.01 and ignore the spin triplet component. The 
subdominant component negligibly changes the spin sus- 
ceptibility, as in the bi- and tri-layer systems. 

In this situation, it is necessary to extend the layer 
dependence of the Rashba spin-orbit coupling, which is 
assumed to be weaker but still existent in inner lay- 
ers. For illustration, we choose two distinct layer depen- 
dences of spin-orbit coupling, as we have no ab initio 
basis for our model. First, we consider a slow reduction 
(screening) of spin-orbit coupling, assuming a decay pro- 
portional to the inverse square of layer distance. Such 
a model yields (ai, a2, as, q;4) = a(l, 9/49, —9/49, —1) 
for the 4-layer system, and (ai, a2, 03, a4, as) = 
a(l, 9/41, 0, -9/41, -1) for the 5-layer system. On the 
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Fig. 11. (Color online) C-axis spin susceptibility in bilayer sys- 
tem for various values while keeping summation -|- = 0.01. 
The ratio changes from to 00. We fix the phase difference 

of the order parameters for the two layers, irrespective of the ra- 
tio |i/)|/|ci|, although interlayer Josephson coupling favors the zero 
phase difference of the dominant order parameter, as in §3. (a) Case 
(A) in which spin singlet order parameter has a 7r-phase difference 
between layers, (b) Case (B) in which spin singlet order parameter 
has a zero-phase difference. We show the results for the interlayer 
coupling = 0.1 and the spin-orbit coupling a = (crosses), 
a = 0.1 (triangles), and a = 0.2 (diamonds). 



ties for different multilayer systems, as shown in Fig. 12. 
We observe a nonmonotonic a-dependence of the spin 
susceptibility because, in our model, spin-orbit coupling 
has several crossover scales (marked by several maxima) 
at ai t±, a2 '-^ t±, and so on. For large a ^ t±, Xs 
recovers the normal state value Xn for even numbers of 
layers, such as 4 and 6, whereas Xs/Xn (M — 1)/M for 
odd numbers of layers, because the center layer has no 
Rashba spin-orbit coupling. 

Next, we consider an extreme situation in which 
Rashba spin-orbit coupling only exists for the outermost 
layers, i.e., a™ = except the outer layers ai = a and 
aM = —a. In this case, the spin susceptibility shows a 
single peak at approximately a ^ t±^ (Fig- 13), as in the 
trilayer systems. For large spin-orbit coupling a ^ t±, 
the spin susceptibility approaches Xs/Xn = 2/M. 

Figures 12 and 13 show nontrivial behaviors of the 
spin susceptibility for multilayer systems, which is af- 
fected qualitatively by the layer dependence of Rashba 
spin-orbit coupling. Recently, the magnetic properties of 
superconducting multilayer cuprates have been investi- 
gated by NMR measurements.^"^' In view of our dis- 
cussion above, it would be interesting to perform similar 
NMR measurements to study the local spin polarization 
and to identify the effect of symmetry-induced spin-orbit 
coupling in multilayer superconductors. 
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Fig. 12. (Color online) Spin susceptibilities along c-axis for M = 
3, 4, 5, 6, and 7. We assume a weak screening effect of Rashba spin- 
orbit coupling, as explained in the text. We assume the uniform 
spin singlet superconducting order parameter ^pm = 0.01 and the 
interlayer coupling tj_ = 0.1. 
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Fig. 13. (Color online) Spin susceptibility along c-axis in more 
than three layer systems with strong screening effect of Rashba 
spin-orbit coupling. Details are explained in the text. The other 
parameters are the same as those in Fig. 12. 



9. Summary and Discussion 

Motivated by the recent investigation of artificially 
fabricated superstructures of CeCoIns/YbCoIns, we 
have theoretically analyzed the basic properties of locally 
non-centrosymmetric superconductors in which spatially 
modulated Rashba spin-orbit coupling plays an impor- 
tant role. The single-particle wave function, supercon- 
ducting gap, and spin susceptibility in the superconduct- 
ing state have been determined. 

Although these layered systems possess a cen- 
ter of inversion symmetry, they display local non- 
centrosymmetricity that affects the spin polarizability 
of the superconducting phase in an interesting way. We 
clearly observe distinct regimes for the response of such 
a layered superconductor with strong or weak interlayer 
coupling. Although the former case exhibits a rather con- 
ventional response in the spin polarization to an exter- 
nal field, the latter refiects properties close to those ex- 
pected for non-centrosymmetric superconductors. An in- 
tuitive understanding of the crossover between the two 
regimes is obtained if we decompose the spin suscepti- 
bility into the conventional Pauli and an additional Van- 
Vleck contributions. The latter results from a snin-orbit 



coupling-induced interband contribution, which is only 
weakly modified by superconductivity, while the Pauli 
part, as an intraband contribution, depends on details 
of the superconducting order parameter. In the discus- 
sion of our model, we demonstrated that the response is 
determined by the interlayer phase structure of the or- 
der parameters, but is independent of the ratio of the 
magnitudes of the singlet and triplet components. 

A special situation appears for dominant spin triplet 
pairing, which has the same phase for all layers. In this 
case, the spin susceptibility along the c-axis is the same 
as that in the normal state and is basically independent 
of spin-orbit coupling. On the other hand, for the dom- 
inant spin singlet pairing, the Pauli contribution to the 
spin susceptibility is completely suppressed at T = 0. 
Here, only the Van-Vleck susceptibility induced by layer- 
modulated Rashba spin-orbit coupling yields a finite con- 
tribution, depending on the interlayer coupling. The de- 
tailed orbital symmetry of order parameters, such as the 
S-, p-, and d-wave, affects these findings only weakly. 

Our study gives the first account of the spin polariz- 
ability of superconductivity in superconductors with lo- 
cal non-centrosymmetricity in clean artificially layered 
superconductors. Previous studies based on similar con- 
cepts have addressed dirty s-wave superconductors^^-* 
and random spin triplet superconductors. ^'^^ Thus, it is 
also interesting to review superconductors with intrinsic 
multilayer structures, such as some high-Tc cuprates, fo- 
cusing on the role of broken local inversion symmetry. 
An additional feature that has been discussed recently is 
the effect of staggered antisymmetric spin-orbit coupling 
on superconductivity in some classes of centrosymmetric 
crystals that also belong to a similar class. 

Finally, we return to the system that initially triggered 
our study, the superlattice of CeCoIn5.^^) In view of the 
fact that bulk CeCoIns is known to realize spin singlet su- 
perconductivity, we believe that the dominantly spin sin- 
glet pairing state is relevant for the multilayer CeCoIus 
in the superlattice. Unfortunately, experimental data of 
spin susceptibility in the superconducting state are not 
yet available. However, the effect of superconductivity on 
the spin susceptibility may be roughly estimated using 
the upper critical field Hc2, because the Hc2 of CeCoIus 
is determined by the paramagnetic depairing effect. From 
a rough estimation, 

Hc2 = , , (52) 

where is the Pauli-limited upper critical field in 
conventional spin singlet superconductors, the large en- 
hancement of the upper critical fields observed in the su- 
perlattice of CeCoIns might be caused by spatially modu- 
lated Rashba spin-orbit coupling and may not necessary 
be a signature of strong-coupling effects. "^^^ Moreover, 
we believe that the variability of the superlattices and 
also the possibility of the local measurements of mag- 
netic properties by NMR measurement would give new 
insights into spin-orbit coupling in these artificial sys- 
tems. Moreover, in this context, the fate of the mag- 
netic quantum critical point of CeCoIns is an issue that 
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arouses experimental and theoretical interest. 
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